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PACS 75.10. Jm - Quantised spin models, including quantum spin frustration 
PACS 75.10.Kt - Quantum spin liquids, valence bond phases and related phenomena 

Abstract - We use the coupled cluster method implemented at high orders of approximation 
to study the spin-1 J 1 -J 2 model on the triangular lattice with Heisenberg interactions between 
nearest-neighbour and next-nearest-neighbour pairs of spins, with coupling strengths Ji > 0 and 
J 2 = reJi > 0, respectively. In the window 0 < re < 1 we find that the 3-sublattice 120° Neel- 
ordered and 2-sublattice 180° stripe-ordered antiferromagnetic states form the stable ground-state 
phases in the regions re < ref = 0.060(10) and re > ref = 0.165(5), respectively. The spin-triplet 
gap is found to vanish over essentially the entire region ref < k < ref of the intermediate phase. 



INTRODUCTION. — Much theoretical progress 
has been made in the search for exotic states in frustrated 
■"O two-dimensional (2D) quantum spin-lattice models that 
Ch have no classical counterparts. A special role is reserved 
0 for the 2D triangular lattice, which is the simplest of the 
1 A rchimedean lattices that exhibit geometric frustration 
( i.e ., which are not bipartite) and the one with the highest 
coordination number. For the triangular-lattice Heisen- 
hercr antiferromagnet (HAFM), comprising spins (with 
spin quantum number s) interacting via a pure nearest- 
neighbour (NN) Heisenberg potential with exchange cou- 
‘Opling Ji > 0, the classical (s —> 00 ) ground state (GS) is a 
3-sublattice Neel state with an angle of 120° between the 
OO spins on different sublattices, thereby breaking the trans- 
CD lational symmetry of the lattice. When Anderson first in- 
troduced the concept of a resonating valence-bond (RVB) 
. .state m, it was hypothesized that the GS of the spin- 
< ^ 1 HAFM on the triangular lattice with NN interactions 
only, could be such an RVB state, the first example of the 
$_i quantum spin-liquid (QSL) states (3J. 

This conjecture was not supported by spin-wave theory 
(SWT), both at lowest (first) and second orders of imple¬ 
mentation, which predict that the system retains the 3- 
sublattice 120° Neel antiferromagnetic (AFM) long-range 
(LRO), albeit with a reduction in the sublattice magneti¬ 
zation (i.e., the magnetic order parameter) of about 50% 
for the s = \ case from the classical (s —»■ 00 ) value Hi]. 
Many other numerical studies, including variational calcu¬ 
lations mm, and some based on the exact diagonaliza- 
tion (ED) of small lattice clusters , lent credence to 

the SWT findings. A variety of recent high-accuracy nu¬ 


merical studies, including some based on the series expan¬ 
sion (SE) method [TS], the ED method m , the Green’s 
function Monte Carlo (GFMC) method [TD], the density- 
matrix renormalization group (DMRG) method |l7j . and 
the coupled cluster method (CCM) [l8liT9l . all now concur 
that the GS of the spin-1 triangular-lattice HAFM with 
NN couplings only has the 3-sublattice 120° Neel ordering, 
with a sublattice magnetization in the rather accurately 
determined range of (40 ± 2)% of the classical value. 

To investigate a possible instability of the 120° Neel 
phase against the formation of more exotic GS phases, it 
is natural to add additional dynamic frustration by, for ex¬ 
ample, including an additional Heisenberg interaction be¬ 
tween next-nearest-neighbour (NNN) pairs of spins with 
an exchange coupling J 2 = re Ji . The study of the resulting 
spin-1 J 1 -J 2 model on the triangular lattice is the main 
aim of the present paper, particularly in the most inter¬ 
esting region of the parameters, J\ > 0 and 0 < re < 1. 

A number of recent studies [TTnnZTl have shown the ex¬ 
istence of a magnetically disordered quantum state in a 
range < n < beyond the range n < for which the 
120° Neel ordering persists, and around the correspond¬ 
ing point = 4 in the classical (s —> 00 ) version of the 
model at which the 120° Neel order melts. Nevertheless, 
there remains disagreement over the values of the quantum 
critical points (QCPs) and and also the nature of 
this intermediate phase. For example, a Schwinger boson 
mean-field theory (SBMFT) study [2D] gives ~ 0.12 and 
0.19. By contrast, two different variational Monte 
Carlo (VMC) studies give k\ 0.05 and ~ 0.18 [21] . 
and ref « 0.10(1) and ref « 0.135(5) [22]. Both VMC 
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Fig. 1: (Colour on-line) The J 1 -J 2 model on the triangular lattice, showing (a) the bonds (Ji = - ; J 2 = -) and the 

Bravais lattice vectors a and 6; (b) the corresponding first Brillouin zone; (c) the 120° Neel antiferromagnetic (AFM) state; (d) 
the infinitely degenerate family of classical 4-sublattice ground states on which the spins on the same sublattice are parallel to 
each other, with the sole constraint that the sum of four spins on different sublattices is zero; and (e) one of the three degenerate 
striped AFM states. For the two states shown the arrows represent the directions of the spins located on lattice sites •. 


studies favour the intermediate state to be a gapless QSL, 
although of different types. However, both the SMBFT 
and VMC techniques involve either uncontrolled approxi¬ 
mations and/or inbuilt bias towards certain forms of order. 

Two essentially unbiased ab initio methods, with con¬ 
trolled approximations, that have been applied to the 
model are the CCM (TS] and the DMRG [231 [251 . Three re¬ 
cent studies give good agreement on the QCPs. The CCM 
study [19] gives rtf ~ 0.060(10) and ss 0.165(5), while 
two separate DMRG studies give k\ « 0.06 and kI] ~ 0.17 
[S2j, and ss 0.08 and n? ~ 0.16 [Mj. (In the third re¬ 
cent DMRG study [25] the model was studied on a three- 
leg cylinder, for which no intermediate phase was found.) 
The two DMRG studies [23][24] both favour the interme¬ 
diate state to be a gapped QSL. Since the nature of the 
intermediate state has not yet been investigated within the 
CCM, which method probably provides the most accurate 
values available of the QCPs kJ and /cfj, our primary aim 
here is to study it via its spin gap. Our main finding will 
be that, in sharp contrast to both recent DMRG studies 
[23I[2I|, the intermediate state is gapless over the entire 
regime < n < refj, just as was found in the two earlier 
VMC studies [2111221 . 

THE MODEL. — We study the Hamiltonian 

H = Ji ^2 Si ■ Sj + J 2 ^2 Si-S k , (1) 

(i,j) «*>k» 

where (i,j) and (( i,k )) run over NN and NNN pairs of 
sites, respectively, counting each bond once only. Each 
site i of the triangular lattice (with lattice spacing a) car¬ 
ries a spin-i particle with spin operator s., = (s?,s/,sf). 
The lattice and the exchange bonds are illustrated in fig. 
Ea). The first Brillouin zone for the triangular lattice is 
a hexagon of side 47r/3a, as shown in fig. E[b). 

Classically, the spin configurations can be described by 
a wave vector Q, such that the spin on site i is given by 

Sj = s[cos(Q ■ rj)ni + sin(Q • rj)n 2 ], (2) 

where ni and h 2 are two orthogonal unit vectors in spin 
space. In the regime J\ > 0, and 0 < k < 1, where k = 


J 2 /J 1 , the corresponding classical (s —> 00 ) phase diagram 
has just one phase transition at kJ 1 = g. For n < the 
system has the 3-sublattice 120° Neel ordering shown in 
fig- Ev), described by either one of the two wave vectors 
representing the two inequivalent corners A and B of the 
first Brillouin zone shown in fig. E,b)> Qa = gp(g,0), 
Q b = For 1 > k > Ki, the classical system 

has an infinitely degenerate family (IDF) of 4-sublattice 
Neel GS phases illustrated in fig. □d), in which the sole 
constraint is + sg + sc + = 0, where Sj denotes 

the spin on each of the four sublattices i = A, B, C , D, 
as shown. The application of lowest-order [be., to 0(l/s)] 
SWT has shown [2M2%] that quantum fluctuations lift this 
accidental degeneracy in favour of the 2-sublattice striped 
states, one of which is shown in fig. Ee), via the well- 
known order by disorder mechanism. These striped states 
have ferromagnetic ordering along one principal direction 
of the triangular lattice and AFM ordering along the other 
two principal directions. They are described by one of 
the three wave vectors representing the three inequivalent 
midpoints K, L, and M of the edges of the hexagonal 
first Brillouin zone shown in fig. [2(b), Q k = ^(0, ^=), 


Q L — a ( 1’ ^3)’ an d Qm — a (■*■> ^3)- 

In our present study we will use both classical GS 
phases, viz., the 3-sublattice 120° Neel phase and the 2- 
sublattice striped state, as reference states, on top of which 
the quantum many-body correlations for the s = g model 
will be systematically incorporated using the CCM. 


THE COUPLED CLUSTER METHOD. - The 

CCM (see, e.g., refs. f29H3l] ) is one °f the most versatile, 
most accurate and most powerful techniques of ab initio 
quantum many-body theory. It has been applied to a huge 
variety of fields in physics and chemistry, including many 
strongly correlated and highly frustrated spin-lattice prob¬ 
lems in quantum magnetism (see, e.g., refs. |18I19I31H37[ ). 
Whereas many alternative methods require a finite-size 
scaling analysis, one of the main strengths of the CCM is 
that it is both size-consistent and size-extensive, and hence 
yields results from the outset in the 2D bulk (N —»• 00 ) 
limit. Another important feature is that it exactly pre- 
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serves the important Hellmann-Feynman theorem at all 
levels of approximate implementation. The method is uti¬ 
lized in practice at various levels of approximation, each 
specified by a truncation index m = 1,2,3,---. The only 
approximation made in using of the CCM is then to ex¬ 
trapolate the sequences of values obtained for physical ob¬ 
servables to the exact m —>• 00 limit. 

Every implementation of the CCM starts with the 
choice of a suitable (normalized) reference (or model) state 
|$), with respect to which the quantum corrections present 
in the exact GS |$) may then be incorporated at the next 
step. Suitable choices for |$) in the present case are the 
two quasiclassical AFM states [viz., the 3-sublattice 120° 
Neel state of fig. [l](c) and the 2-sublattice striped state of 
fig.m e)]. The exact ket and bra GS wave functions |\I/) 
and (*1/ |, which solve the respective Schrodinger equations 
H |*F) = £7|\I/) and ($|if = i£($|, are chosen to have nor¬ 
malizations such that (’Flu/) = (<I>|\I/) = ($|$) = 1. In 
terms of the chosen model state the CCM now utilizes its 
distinct exponential parametrizations, |\H) = e s |$) and 
(*j/| = (<F| ( Se~ s '. The two GS correlation operators are 
themselves formally decomposed as S = J2i^oSiCf and 
S = 1 + SiCj, where we define C$ = 1 to be the 

identity operator, and where the set index I denotes a 
complete set of single-particle configurations for all of the 
N particles (be., in our case, spins). 

For spin-lattice applications each set index I represents 
a unique multispin-flip configuration with respect to |$), 
such that the corresponding wave function for this configu¬ 
ration of spins is Cf |$). We thus require |$) to be chosen 
as a fiducial (or cyclic) vector (be., in physical terms, a 
generalized vacuum state) with respect to the complete 
set of mutually commuting many-body creation operators 
{Cj~}. In other words, we require Cj~ and its destruc¬ 
tion counterpart, CJ = (C/)I, to satisfy the conditions 
($|C+ = 0 = C7|$), V/^0. 

It is useful to treat each lattice site i as equivalent 
to all others, whatever the choice of |$). To that end 
we passively rotate each spin Sj so that in its own lo¬ 
cal spin-coordinate frame it points in the downward (be., 
negative z s ) direction. Henceforth our description of 
the spins is given in terms of these locally defined spin- 
coordinate frames, in which all independent-spin product 
model states take the universal form |$) = | 4-4-4- • ■ ■ 4- 
). Similarly, Cf may be written in the universal form 
Cf = sjs^ • • • sjj - ; n = 1,2, ■ ■ • , 2 sN, as a product of 
single-spin raising operators, s k = s k +is v k . The set index 
/ = {hi, & 2 , ■ ■ • k n \ n = 1,2, ■ • • , 2 sN} thus becomes a set 
of (possibly repeated) lattice site indices. For the present 
case, with s = each site index included in set index / 
may appear no more than once. Once the local spin coor¬ 
dinates have been chosen for a given |$), the Hamiltonian 
H needs to be re-expressed in terms of them. 

The GS CCM c-number correlation coefficients {Si, Si} 
are now calculated by requiring that the GS energy ex¬ 
pectation functional, H = H[Si,Si] = ($|S'e _,s ibe ,s |$), 


be minimized with respect to each of them sepa¬ 
rately. This yields the coupled set of nonlinear equa¬ 
tions {$\CJe~ s He s \<&) = 0; V/ 0, for the coefficients 
{Si}, and the coupled set of linear generalized eigenvalue 
equations ($|S'[e _ ' s ibe' s , C/]| < I>) = 0 ; or equivalently, 
($| S{e~ s He s — E)Cf |$) = 0, VJ ^ 0, for the coefficients 
{Si}, once the coefficients {Si} are solved for. Having 
solved for {Si, Si}, one may then, for example, calculate 
all GS properties such as the energy E = ($|e^ s ibe s |$) 
and the magnetic order parameter (be., the average local 
on-site magnetization) M = — jf(®\SJ2k=i e _,s s^e ,s | < I)). 

Excited-state (ES) wave functions are parametrized 
within the CCM as |$ e ) = X e e s \<&), in terms of a lin¬ 
ear excitation operator X e = J2i^o^i^t ■ By combin¬ 
ing the GS Schrodinger equation with its ES counterpart, 
H |$ e ) = E e |$ e ), and using that [ X e ,S} = 0, we find 
e~ s [H, X e ]e s |$) = A e A e |$), where A e = (E e — E) is the 
excitation energy. By taking the overlap with the state 
($|CJ we find {<f>\Cj [e~ s He s , X e }\<!>) = A e A e |$); V/ ^ 
0, where we have used that the states labelled by the 
set indices / are orthonormalized, |$) = 8{I, J ). 

These generalized eigenvalue equations may then be solved 
for the set {Xf} and A. 

No approximations have yet been made. One might ex¬ 
pect that truncations are needed in the evaluation of the 
exponential functions e ±s in the CCM parametrizations. 
However, we note that these always appear only in the 
form of the similarity transform e~ s He s of the Hamilto¬ 
nian in all of the GS or ES equations we need to solve. 
This similarity transform may be expanded in terms of 
the well-known nested commutator series. Another key 
feature of the CCM is that this otherwise infinite series 
actually now terminates exactly at the double commuta¬ 
tor term. This is due to the basic SU(2) commutation 
relations and the fact that all of the terms in the expan¬ 
sion EimSiC+ for S commute with one another and 
are simple products of single-spin raising operators. Fur¬ 
thermore, all terms in the expansion for e~ s He s are thus 
linked. Hence, the Goldstone linked cluster theorem and 
the consequent size-extensivity of the method are exactly 
preserved even if the expansion for S is truncated. 

Thus, the only approximation made in practical imple¬ 
mentations of the CCM is to restrict the set of indices {/} 
retained in the corresponding expansions for the operators 
{S', S'}. As in our previous work on this model [T5], and 
in many other applications of the CCM to quantum mag¬ 
nets, we use here the well-tested localized (lattice-animal- 
based subsystem) LSUBm scheme. At the mth level of 
approximation all multispin-flip configurations {/} are re¬ 
tained in the expansions for S and S that are constrained 
to no more than m contiguous sites. A cluster configu¬ 
ration is contiguous if every site in it is NN to at least 
one other. By utilizing the space- and point-group sym¬ 
metries of the model, together with any conservation laws 
that apply to both the Hamiltonian and the model state 
under study, we may reduce the terms retained to a min- 
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Fig. 2: (Colour on-line) CCM results for the GS magnetic order parameter M versus the frustration parameter k = J 2 /J 1 , for 
the spin-4 J 1 -J 2 model on the triangular lattice (with J\ >0). The left and right sets of curves in each panel are based on the 
120° Neel and striped AFM states, respectively, as CCM model states. The LSUBm results are shown: (a) with m = {3, 5, 7, 9}; 
and (b) with m = {4, 6, 8,10}. In both cases the extrapolated LSUBoo results using the respective data set in eq. © are shown. 
We also show with cross (x) symbols the corresponding extrapolated values based on the 120° Neel state using eq. © for the 
case k = 0 (i.e., the pure triangular-lattice HAFM with NN interactions only). 


imal number Nf = Nf(m) of fundamental LSUBm con¬ 
figurations. Nevertheless, Nf(m ) increases rapidly as the 
truncation index m is increased, and one soon needs to 
use massive parallelization and supercomputing resources 
for the higher-order calculations [321138] . 

For the ES calculation of the (magnon, triplet) spin 
gap A the choice of cluster configuration {/} retained in 
the expansion for the excitation operator X e is different 
from those for the GS correlation operators S and S, and 
different too for each model state, since we now restrict 
attention to those that change the z component of total 
spin, S z , by one unit (in the local spin-coordinate frames), 
rather than to those with S z = 0 for the GS calculation. 
However, to ensure comparable accuracy for both the GS 
and ES calculations, we use the LSUBm scheme in each 
case. The number Nf(m), for a given value of m, is appre¬ 
ciably higher for the ES case than for the corresponding 
GS case using the same model state. The corresponding 
GS CCM LSUB?n equations have been solved for both 
model states for m < 10. Whereas we can also solve the 
ES LSUBm equations for m < 10 for the striped state, 
the number Nf( 10) for the ES using the 120° Neel state 
is now sufficiently large that we are restricted in practice 
to m < 9 in this case. 

Clearly, the CCM LSUBm approximations become ex¬ 
act, by construction, in the m —> 00 limit. Thus, although 
the CCM obviates the need for any finite-size scaling, we 
do need to extrapolate our LSUBm results to the m —> 00 
limit as our last step and as our sole approximation. There 
exist very well-tested accurate extrapolation rules for such 
GS quantities as the energy per spin, E/N , and the mag¬ 
netic order parameter, M, as we have described and used 
in our previous paper m on this model. For example, as 
we discussed there, for unfrustrated models or ones with 


relatively minor amounts of frustration we use 

M(m) = b 0 + b\m~ l + 6 2 m -2 , (3) 

whereas for highly frustrated systems, especially ones close 
to a QCP, we use 

M(m) = bo + frimT 1 / 2 + 62 m -3 / 2 . (4) 

For the spin gap we use the extrapolation scheme mm 

A(m) = do + dim ” 1 + ^ui ” 2 . (5) 

Of course, for the GS expectation value of any parameter 
P, we may always fit our LSUBm data to an asymptotic 
form P(m) = po + piin~ v , in order to extract the leading 
exponent v. For our A(m) results, for example, we find 
v ss 1 for all values of k, which hence justifies the use of 
Eq. ©. Whenever possible we fit each of the extrapo¬ 
lation schemes of eqs. © ©, all of which contain three 
fitting parameters, to LSUBm results with at least four 
different values of m. Since the triangle is the basic struc¬ 
tural element of the lattice, the lowest-order results with 
m < 3 are excluded from the extrapolations. 

RESULTS. — We show in fig. [2] our results for the 
magnetic order parameter M for the model, using both 
the 120° Neel and striped AFM states as CCM model 
states. Since the LSUBm results show a marked even-odd 
staggering effect in the truncation index m, we display the 
odd and even values separately in figs, ©a) and[ 2 ](b), to¬ 
gether with the associated LSUBoo extrapolations using 
eq. @, which is appropriate in the highly frustrated re¬ 
gions near the QCPs, together with the value using eq. 
© for the unfrustrated case k = 0 only. The QCPs (viz., 
the points n at which the extrapolated LSUBoo values for 
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Fig. 3: (Colour on-line) CCM results for the spin gap A versus the frustration parameter k = J 2 /J 1 , for the spin-1 J 1 -J 2 
model on the triangular lattice (with Ji = 1). The left and right sets of curves in each panel are based on the 120° Neel and 
striped AFM states, respectively, as CCM model states. The LSUBm results are shown: (a) with m = {3, 5, 7, 9}; and (b) with 
m = {4, 6, 8} for the 120° Neel state and m = {4, 6, 8,10} for the striped state. In both cases the extrapolated LSUBoo results 
using the respective data sets in eq. © are shown, together with error bars associated with the fits using four data points. 


M vanish) show great consistency between those based on 
even m values and those based on odd m values. A care¬ 
ful analysis of our errors gives our best estimates for the 
two QCPs as k\ = 0.060 ± 0.010 and k c 2 = 0.165 ± 0.005. 
We note, parenthetically, that the smooth shape of the 
M{k) curves in fig. [2] near both QCPs lends support to 
both transitions being continuous. However, since the ex¬ 
trapolations near any critical point are very sensitive, we 
cannot entirely rule out the possibility of either transition 
being first-order. 

An interesting feature of our LSUBm results displayed 
in fig. [2] is that we find solutions based on a given model 
state over a region that extends beyond the associated ex¬ 
act (LSUBoo) QCP, typically out to some associated ter¬ 
mination point K*(m), beyond which, at a given LSUBm 
level of approximation, no real solution exists. Such ter¬ 
mination points have been found in many other applica¬ 
tions and are very well understood. They have been shown 
[31] to be direct manifestations of the associated QCP (at 
k c , say) at which the order of the model state melts in 
the physical system. Typically the extent | K 4 (m) — k c \ of 
the region of unphysical intrusion of the solution with the 
order properties of the CCM model state into the neigh¬ 
bouring phase (beyond the QCP at n c ) in which this form 
of order vanishes, decreases as the truncation index m is 
increased, and vanishes in the exact limit m —> oo. It is 
clear from fig. [2] that these regions of unphysical intrusion 
of the CCM LSUBm solutions into the regime < k < k 2 
of the intermediate phase cover essentially the whole in¬ 
termediate region for all values m < 10 on both the neigh¬ 
bouring 120° Neel and striped quasiclassical AFM sides. 
For that reason we are now well placed to investigate the 
intermediate phase regime using both quasiclassical model 
states. 

Thus, in figs. Ha) and ©b) we show our correspond¬ 


ing LSUBm results for the spin gap A, together with the 
extrapolated LSUBoo values using eq. ©. Goldstone’s 
theorem implies that any state breaking spin-rotational 
symmetry, including all states with magnetic LRO, must 
have a vanishing gap. We see from figs, ©a) and©b) that 
our extrapolated LSUBoo results for A are wholly consis¬ 
tent with A = 0 for the striped state, k > k 2 « 0.165. 
The results on the (computationally more challenging) 
120° Neel-ordered side, k < « 0.06, are slightly less 

accurate, but still compatible with Goldstone’s theorem, 
within small errors, presumably associated with extrapo¬ 
lations in this region being less accurate. 


Most interestingly, the results in fig. [3] based on both 
quasiclassical AFM states show no tendency at all for the 
spin gap A to increase in the regions of the intermediate 
phase accessible by the unphysical intrusion effect. Of 
course, the spin gap, by definition, satisfies A > 0. The 
regions where A < 0 in fig. [3] are clearly wholly unphysical 
therefore. However, our results are invalid in the regions 
beyond where the curves based on each model state cross, 
since the value of A is, by definition, unique. Beyond the 
crossing point we are simply too far into the unphysical 
regime, and we see from fig. [3] that this is reflected in the 
large error bars there. We note too that our largest error 
for A on the 120° Neel ordered side (where A = 0 is exact) 
is A < 0.1, which is compatible with the values we obtain 
for A in the intermediate regime. The extrapolated values 
of the spin gap are thus compatible with the spin gap 
vanishing essentially everywhere within the intermediate 
phase region, < k < k 2 . For comparison, we note that 
the estimate obtained by a recent DMRG study [23] for 
the bulk triplet gap is A ss 0.3 ± 0.1 at k = 0.1, which is 
clearly at odds with our results. 
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SUMMARY AND DISCUSSION. - Our primary 
finding is that we find no evidence at all for any part of the 
intermediate phase to be a gapped state. We may com¬ 
pare our study with similar recent CCM spin-gap studies 
of both the spin-^ J1-J2 model on the square lattice [591 
and the spin-^ J\ J2-J3 model on the honeycomb lattice 
(with J 3 = J 2 ) [40 • For the former model A was again 
found to be zero over that part of the equivalent inter¬ 
mediate phase region accessible by the unphysical intru¬ 
sion mechanism, which was not the whole region in this 
case. By contrast, for the latter model, strong evidence 
was found that A / 0 over the similarly accessible part 
of its equivalent intermediate phase, which finding pro¬ 
vides strong support for that phase having valence-bond 
crystalline (VBC) order, as is now broadly agreed. 

Finally, in contrast to the recent DMRG studies mm 
that find the intermediate phase of the spin- 4 J\-J 2 model 
on the triangular lattice to be a gapped QSL, our CCM 
results indicate that the intermediate phase is gapless. 
While our results clearly preclude the intermediate phase 
in this model to have any form of VBC order, since any 
such state is gapped, they point towards a gapless QSL as 
the most likely candidate. 
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